Numerical Construction of LISS Lyapunov Functions under a 

Small Gain Condition 

CN Roman Geiselhart Fabian Wirth* 

^ January 20, 2013 

00 



< 



Abstract 



In the stability analysis of large-scale interconnected systems it is frequently desirable 
to be able to determine a decay point of the gain operator, i.e., a point whose image 
under the monotone operator is strictly smaller than the point itself. The set of such 
decay points plays a crucial role in checking, in a semi-global fashion, the local input-to- 
4^ state stability of an interconnected system and in the numerical construction of a LISS 

Lyapunov function. We provide a homotopy algorithm that computes a decay point of a 
monotone operator. For this purpose we use a fixed point algorithm and provide a function 
whose fixed points correspond to decay points of the monotone operator. The advantage 
to an earlier algorithm is demonstrated. Furthermore an example is given which shows 
how to analyze a given perturbed interconnected system. 

Keywords: homotopy algorithm, monotone operator, LISS Lyapunov function, interconnected sys- 

tern, small gain condition 

In recent years large-scale systems have received renewed attention with applications in for- 
^ mation control, logistics, consensus dynamics, networked control systems and further appli- 

1-H cations. While stability conditions for such large-scale systems have already been studied in 

[20, 26, 32] based on linear gains and Lyapunov techniques, nonlinear approaches are more 
^ recent. The groundbreaking concept that has proven fruitful is the notion of input-to-state 

^ stability (ISS) as introduced in [28]. 

^ For large-scale nonlinear systems it may be difficult to prove ISS directly, but if a large-scale 

system is defined through the interconnection of a number of smaller components, which 
are ISS, then there exist small gain type conditions guaranteeing the ISS property for the 
interconnected system. For the case of two subsystems this result was obtained in [15, 14] 
both in a trajectory based as well as a Lyapunov formulation. Recently, there has been a 
substantial effort to extend these results to the case of a greater number of subsystems, see 
[7, 8, 12, 18, 22, 4]. It is the purpose of this paper to provide numerical methods that make 
some of the available results applicable for practical problems. 

The general setting is here to consider a number of systems that are input-to-state stable with 
respect to external and internal inputs. The effect of the subsystems, described by comparison 
functions, is collected in the gain matrix T. The special structure of the interconnected 
system now leads to a monotone operator T^ on the positive orthant W^. So-called monotone 
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aggregation functions can be used to formulate the effect of several inputs on a system in a 
general manner. Standard examples of such functions are summation and maximization, but 
in [8] some examples are provided that also other types of aggregation functions may be useful 
depending on the system under consideration. We would like to point out that the particular 
relation of the maximization and summation formulation of small gain conditions is analyzed 
in [5]. In [16] the authors study interconnections where small gain conditions are satisfied 
after certain transient periods and derive stability results. 

Many available small gain results state that input-to-state stability for the overall system 
follows from the existence of a so-called 0-path with respect to T^, [4, 5, 6, 7, 8, 21]. Further- 
more an ISS Lyapunov function for the interconnected system can be constructed using this 
path and the ISS Lyapunov functions of the subsystems. Note that also for other small gain 
type formulations as the cycle condition in the maximization case or the spectral radius con- 
dition in the linear summation case, it may be seen that these conditions can be equivalently 
formulated in terms of Jl-paths. 

In [8] the construction of an fi-path is described. The crucial ingredient that usually cannot 
be obtained in a straightforward manner is a decay point of F^, that is a point s G M+ for 
which F^(s) <^ s in the order induced by the cone W^. Once such a point is found there 
are straightforward numerical procedures for the construction of Lyapunov functions or for 
checking the ISS property. There are two particular cases in which a straightforward way 
is known to compute decay points: (i) If the gains are linear and summation is used, then 
the problem becomes one of checking whether the spectral radius of F is below 1 and finding 
an appropriate eigenvector. As F is nonnegative this problem is particularly easy and well 
studied; (ii) If the maximization formulation of ISS is used then a very nice observation of 
[17] is that, provided a small gain-condition holds, for s* := max{s, F(s), . . . , Fjy^^(s)}, we 
have F(s*) < s* for all s G M.^, which is almost a decay point. The methods presented in this 
paper are suitable for the cases that the problem at hand is not within one of the two classes 
described above. In this paper we provide numerical procedures for computing such points 
and thus also for local O-paths. We call the approach semi-global because it does not require a 
priori restrictions. In particular, if a small gain condition is satisfied globally, then the design 
variables of the algorithm can in principle be chosen so that the numerically guaranteed region 
of stability is arbitrarily large. 

As we compute a decay point numerically the overall construction of Lyapunov functions as 
well as the verification of the ISS property is only performed locally. Indeed, the approach 
relies on local results of the small gain type. Local small gain results have been considered 
in [3, 2, 13] in an input-output operator context, resp. for discrete-time systems. In [6] local 
ISS (LISS) definitions and local small gain theorems within the framework considered here. 
In this work the knowledge of a decay point leads to the local input-to-state stability of the 
interconnected system and to the construction of a LISS Lyapunov function. 
The algorithm developed here, that computes a decay point for a given monotone operator F^, 
is a particular simplicial fixed point algorithm (SFP-algorithm) customized in such a way that 
we obtain a decay point of F^. To ensure the convergence of the SFP-algorithm we require 
irreducibility of the gain matrix F. This is no significant restriction because by standard graph 
theoretic algorithms the irreducible components of the system can be obtained efficiently, [30] . 
The paper is organized as follows. In Section 1 we provide the necessary notions and a 
short introduction to comparison functions and graphs. In Section 2 we recall the Lyapunov 
formulation of ISS for interconnected systems, give a local small gain theorem and outline the 
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construction of a LISS Lyapunov function for the overall systems. Section 3 contains the main 
results of this paper. First we recall some facts about homotopy algorithms and introduce 
the SFP-algorithm where we mainly follow the book of [33]. In subsection 3.4 we state some 
sufficient conditions on and prove that the SFP-algorithm converges to a decay point of F^. 
At the end of this section some improvements of the algorithm are discussed. We conclude 
this work in Section 4 where we discuss two examples. The first one shows that this new 
algorithm improves on an earlier algorithm that is due to a homotopy algorithm of Eaves [9] 
(cf. [24]) where we revisit a nonlinear example from [23]. In the second example we use our 
techniques to show numerically that a particular perturbed interconnected system is LISS. 



1 Preliminaries 

1.1 Notation and conventions 

Let M denote the field of real numbers, M-|- the set of nonnegative real numbers, and M.^ (resp. 
M^) the vector space of (nonnegative) real column vectors of length N. Then induces a 
partial order for vectors v,w E M^. We denote v > w <^=^ Vi > Wi, v > w <^=^ v > w and 

V ^ w, V ^ w <^=^ Vi > Wi, each for i = 1, . . . ,N, where Vi denotes the i*^ component of 
the vector v. Let v,w ^ M.^ be given. Then we define the order intervals [v, w] := {x G M.^ : 

V < X < w} a V < w, {v, vj) := {x G : v <^ x <^ w} ii v <^ w, and analogously the order 

intervals {v,vj] and [v,w). For x G M.^ we use the Euclidean norm ||a;|| = \Jj2^i kiP- The 
space of measurable and essentially bounded functions is denoted by L°° = L°°([0, oo);M^) 
with norm 11 • ||oo. 



1.2 Comparison functions and induced monotone operators 

To state the stability definitions that we are interested in, three sets of comparison functions 
are used. We call a function a : M+ — >■ M_|_ a function of class /C, if it is strictly increasing, 
continuous, and satisfies a(0) = 0. If a G /C is unbounded, it is said to be of class /Coo- A 
function /3 : M_|_ x M_|_ — >■ M+ is called a function of class KC, if it is of class /Coo in the first 
argument and strictly decreasing to zero in the second argument. It is easy to see that if 
p G /Coo, then its inverse p^^ : M+ — )■ M+ exists and is also of class /Coo. 
To formulate general small gain conditions we need the following definition, see [8]. 

Definition 1.1 A continuous function p : M.^ — >■ M_|_ is called a monotone aggregation func- 
tion, if the following properties hold: 

(i) positivity: p{s) > for all s G and p{s) = 0, if and only if s = 0; 

(a) strict increase: p{s) < p{t), if s <^ t; 

(Hi) unboundedness: p{s) — >■ oo, if \\s\\ — >■ oo. 

The space of monotone aggregation functions is denoted by MAF^r. 

The properties in Definition 1.1 can be extended to vectors in the sense that p = {pi, . . . ,p]\[)~^ G 
MAF^, Pi G MAFN,i = l,...,N, defines a mapping from M^^^ to by p{A)i = 
Piian, Oijv) for A = {aij)f^j^^ G R^""^ . 
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We want to generalize this to matrices of the form T = (7ij)i^=i G (/Cqo U {0})^^^, where 
denotes the zero function. This leads to an operator : defined by 

/ Ati(7ii(si)> • • • >7iiv(sAr)) \ 
r^(s) := in o r)(s) := : G for s G R^. (1) 

\ Mw(7Ari(si), . . . ,7Ariv(sAr)) / 

For the A; times composition of this operator we write F^. We call F^ 

(i) monotone, if F^('u) < Tij,{w) for all u; G R^ with v < w; 

(ii) strictly increasing, if F^(f) ^ r^(tt;) for all v,w E R^ with v <^ w. 

Remark 1.2 iVote t/iat if T e (/Coo U {0})^^^ and /x G MAFj^, t/ien F^ is monotone and 
satisfies F^(0) = 0. 

The next definition is fundamental in the following. 

Definition 1.3 For a given function T : M.^ — t- we define the set of decay by 

n{T) := (s G : T{s) < s} . 

For short we just write i7, if the reference to T is clear from the context. Points in Q are 
called decay points. 

1.3 Graphs and matrices 

A directed graph G{V, E) consists of a finite set of vertices V and a set of edges E <ZV xV . 
If G{V, E) consists of N vertices, then we may identify V = {1, . . . , N}. So if (j, i) G E, then 
there is an edge from j to i. The adjacency matrix Aq = (oij) of this graph is defined by 
Oij = 1, if {j, i) & E and aij = else. We call the graph G{V, E) strongly connected, if for each 
pair (i, j) there exists a path {eiQ,n,ei^^i^, . . . , ei^^_^^i^) with i = io,j = ik such that ej(_-^,i( G E 
for all i = 1, . . . , A;. It is well known that the graph Giy, E) is strongly connected, if and only 
if the adjacency matrix Aq is irreducible, i.e., there exists no permutation matrix P such that 

A-P^(^ ^ 

for suitable, square matrices B and D. These definitions can be carried over to matrices 
F G (/Coo U {0})'^^'^. To this end we define the matrix Ay = (aij)fj^i by a^j = 1, if jij G /Coo, 
and Qij = 0, if jij = 0. We call F irreducible, if the matrix is. 

2 Input-to-state stability and small gain theorems 

Consider the control system 

x{t) = f{x{t),u{t)), teR+, (2) 

where u G R™ is the input and x G M" is the state. We assume that / : R" x R" is 

continuous and locally Lipschitz in x uniformly for u in compacts; by this we mean that for 
every compact Ki C M" and compact subset K2 C R"^ there is some constant c > such that 
||/(x, u) — f{z, u)\\ < c\\x — z\\ for all x,z e Ki and all u G K2. Further we assume /(0, 0) = 
and all solutions can be extended to [0,oo). 
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Definition 2.1 Consider the system (2) and let V : — )■ M+ 6e continuous and locally 
Lipschitz continuous on M"\{0}. Then V is called an ISS Lyapunov function, if there exist 
ai,a2 G A^oo suc/i that for all x G M", 

ai(|N|)<F(x)<a2(|N|), (3) 

ancZ «/ there exist 7 G /C ancZ a positive definite function such that for all u G and almost 
all X e W, 

V{x) > j{\\u\\) VV{x)f{x,u) < -«3(lk||). (4) 

Note, that we only assume Lipschitz continuity of the ISS Lyapunov function V. By Rade- 
macher's Theorem, see e.g. [10], this imphcs that V is differentiable almost everywhere and 
we consider the decay condition (4) only at points where V is differentiable. An equivalent 
formulation can be given in terms of Clarke subdifferentials but we refrain from doing so, since 
this will play no further role in the paper, sec also [7, 8]. 

System (2) is called input-to-state stable (ISS), if it has an ISS Lyapunov function. There 
is another, trajectory-based definition of ISS which is equivalent to the existence of an ISS 
Lyapunov function (cf. [29] for smooth Lyapunov functions and [8, Theorem 2.3] for continuous 

and locally Lipschitz continuous functions). 

Now we want to generalize this stability definition to networks. Let N and consider the 
A'^ interconnected systems given by 

Xl = fi{xi,...,XN,u) 

\ ■ (5) 

XN = fN{xi,. . ■ ,XN,u) 

Assume that Xi G R"-', u G M"^ and the functions fi : R^j=i — > arc continuous and 
locally Lipschitz in x = {xj , . . . ,xjf)^ uniformly for u in compacts. Let Xi denote the state 
of the i^^ subsystem and assume « as an external control variable. Without loss of generality 
we may assume to have the same input for all systems, since we may consider u as partitioned 
u = {uj , . . . , ujf)'^, such that each Ui is the input for subsystem i only. Then each fi is of the 
form fi{. . . ,u) = fi{. . . , TTi^u)) = fi{. ..,Ui) with a projection tTj. 

If we consider individual systems, we treat the state Xj,j / i, as an independent input for Xj. 
Assume that for each subsystem i e {1,. .. , N} there exists a continuous and locally Lipschitz 
continuous function Vi : M"' — )• M+ such that for suitable an, a2i G /Cqo 

aiii\\xi\\)<Viixi)<a2ii\\xi\\) for all G M"' . (6) 

We call the function Vi an ISS Lyapunov function for the subsystem i, if there exist iii G 
MAFjv+i, Jij G /Coo U {0},j 7^ i,7j„ G /C U {0} and a positive definite function ccj such that 

Vi{Xi) > fli{jil{Vi{xi)),...,'yiN{VN{xN)),liu{\\u\\)) 

=> VVi{xi)fi{x,u) < -aiiWxiW). (7) 

The functions jij and 'fiu are called ISS Lyapunov gains. We distinguish between the internal 
inputs Xj and the external input u of the i*^ subsystem. These gains indicate the influence 
of the inputs on the state. This is why we set 7^^ = 0, if fi does not depend on Xj and we 
collect the internal inputs into the gain matrix T := (7jj)^-^^. Note that V and the /Xj define 
a monotone operator : as in (1) (cf. Remark 1.2). 
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2.1 A local small gain theorem 

In this section we assume that the interconnected system (5) satisfies an ISS condition of the 
form (7) for ISS Lyapunov functions V^, i = 1, . . . , iV. Denote the corresponding gain operator 
by as in (1). We assume that F is irreducible, so that is strictly increasing (cf. [21, 
Lemma 2.7]). A local ISS Lyapunov function for the overall system given by 

x = f{x,u) (8) 

and X = (xj , . .., x]^)'^ , f = {fj, . .. , fjf)'^ may now be constructed as follows. 
Assume there exists a u; S> with 

r^{w) < w. (9) 

Then the sequence r^{w), A; = 1, 2, ... is strictly decreasing and so limfe_>oo ^^{w) exists. If 



lim T'^iw) = 0, (10) 



then we define the linear interpolation of the points {r^(T^)}fcgN by a : [0, 1] M^: 

Jo, ifr = 

^^^^ = + k) [ll - r\V%w) + [r - ^]r^-n^^)) , if r G \l k e N. ^^^^ 

Note that a is continuous on [0, 1] by (10) and strictly increasing in all component functions 
as is assumed to be irreducible. With this construction local ISS Lyapunov functions can 
be constructed using the following summary of existing results (cf. [6, Theorem 5.5]). 

Theorem 2.2 Assume that system (5) satisfies ISS conditions of the form (7) for all i = 
1,. . . ,N, and that the gain matrix T is irreducible. If there exists an w ^ so that (9) and 
(10) hold, then a local ISS Lyapunov function for the overall system (8) is given by 

V{x)= max a;\Vi{xi)). (12) 
i=l,...,N ^ ^ ^ " 

In particular, the implication 

V{x)>^{\\u\\) => VV{x)-f{x,u)<-a{V{x)) (13) 
holds locally with 7 G /Cqo given by [6, Proposition 4-3]. 

Rem£irk 2.3 (i) By "local" we mean "in an open neighborhood of the origin (x*,u*) = 
(0,0)". In particular, [6, Theorem 5.5] shows that the assured domain of stability in- 
creases with the choice of w. In particular, if the small gain condition holds globally the 
domain where (13) holds can be made arbitrarily large. 

(a) By (9) and (10) we have the small gain condition r^(s) ^ s for all s G [0, f/;]. 

(Hi) Note that a{r) G J^(ry^) for all r G [0,1] and a belongs to the class of ^ -paths (cf. [8, 
Definition 5.1]). 
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Remark 2.4 (i) Theorem 2.2 is the starting point for our numerical considerations. If we find 
the decay point w, then the problem of constructing Lyapunov functions or checking small gain 
conditions becomes easy. In the remainder of the paper we concentrate on giving numerically 
tractable solutions to this problem. 

(a) In the linear case with fi = T, we have rs(s) = Ts with T G M^^^. Here the existence 
of a decay point w ^ with Tw <^ w is equivalent to the spectral radius of T being less than 
one, i.e., 1 > p(r) = {|A| : X is an eigenvalue ofT} (cf. [21, Lemma 1.1]). So finding a decay 
point is just an eigenvalue problem. This is why we assume F^^ to be nonlinear. 

3 A homotopy algorithm for computing a decay point w G ^^(r^) 

In this section we want to develop an algorithm that computes a decay point w G ^(Tf^) for 
a given continuous and monotone operator : M^^. We know that such a point exists 

for any norm, if the small gain condition 

r^{s) t s for all s G M^\{0} (14) 

is satisfied (cf. [7, Proposition 5.3]). 

To find such a point we will extend a homotopy algorithm that was also used by Merrill [19] to 
compute fixed points of upper-semicontinuous (u.s.c.) point-to-set mappings. Note that since 
a continuous single- valued function is in particular an u.s.c. point-to-set mapping our problem 
falls in the class of problems that can be treated by homotopy algorithms. However, Merrill's 
condition introduced in [19] is not sufficient for convergence in our case, as the domain of the 
mapping is only the nonnegative orthant. The idea to the design of a convergent algorithm is 
to construct a function ^ : M.^ — >■ M.^ , which has the property that its fixed points are decay 
points of r^, and to show that the homotopy algorithm will converge to approximate fixed 
points of (p, which arc also decay points of T^. This algorithm is semi-global since by choosing 
design variables appropriately we end up in a decay point with arbitrarily large norm. 
In Section 3.1 we present the triangulation we need for the computation of fixed points. Before 
introducing the homotopy algorithm in Section 3.3 we first provide some facts about homotopy 
algorithms in Section 3.2. In Section 3.3 we mainly follow the book of Yang [33, Section 4.3]. 
In Section 3.4 we will give the function mentioned above and prove the convergence of the 
SFP-algorithm. Finally in Section 3.5 we give approaches for improving the algorithm and 
further give suggestions for the choice of the design variables used in the mapping 0. 

3.1 Simplices and triangulations 

We briefiy recall facts about covering convex sets by triangulations. A set C C M.^ is called 
convex, if for all a, G C it holds Sa,b '■= {Aa + (1 — A)6 : A G [0, 1]} C C. The convex hull 
co{M} of a set M C is the smallest convex set containing M. If M is finite, we also say 
that co{M} is spanned by M and denote this by ({f* G M}). The dimension of a convex set 
is equal to the dimension of the smallest affine subspace U C containing C. 

Definition 3.1 An -simplex S is an N -dimensional, convex polytope spanned by N + 1 
vectors v^,.. .,v^+^ in M^, M >N, i.e.. 
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A subsimplcx q of S is a simplex spanned by a subset of the set of vertices of S, i.e., <^ = 
{w* : i G /c;} with /,j C {1, . . . , + 1}. Zero -dimensional subsimplices are just the vertices of 
the simplex, one- dimensional subsimplices are called edges between the vertices and {N — 1)- 
subsimplices are called facets. The subsimplex S{j) = ({f* : i / j}) is called the facet opposite 

Clearly, since any A^-simplex is A^^-dimensional, of the vertices are linearly independent 
and it holds ^ for i ^ j. Simplices can be used to cover convex sets in as follows. 

Definition 3.2 Let C be an m-dimensional convex set in M^. A set T of m-simplices is 
called a triangulation of C, if 

(i) C is the union of all simplices in T; 

(ii) for any ryi, r/2 G T, r/i ^ ri2, the intersection 771 fl 772 is either the empty set or a common 
facet of both; 

(iii) every x E C has an open neighborhood intersecting only a finite number of G T. 

By T'^ we denote the set of all A;-subsimplices of T. It is easy to see that = T and 7^ de- 
scribes the set of the vertices of the simplices in T. To distinguish simplices, or triangulations, 
we introduce the diameter of a simplex G T by 

diam(77) = max{||x — y|| : x, y G rj} 

and the mesh size of a triangulation T by 

mesh(T) = sup{diam(?7) : ry G T}. 

There is one special triangulation of M^, which will be used to compute decay points. Let 
denote the z*'* unit vector in R^. The iCi -triangulation is defined as the set of all A?^-simplices 
with vertices x^, . . . , x^~^^ such that 

x^ G and x'+^ = x' + e^^(j) for alH G {1, . . . , N}, 

where ttn = (7rAr(l), . . . , 7riv(-/V)) is a permutation of the elements of the set {1,...,A^}. 
We denote these simplices by 77(3;^, ttat)- See [33, Theorem 1.4.8] for a proof that Ki is a 
triangulation in the sense of Definition 3.2. An illustration of this triangulation is given in 
Figure 1. 

Defining SC = {6x : x e C} hi C C , 6 > 0, and 6F = {5C : C e F} for a family F of 
subsets of we obtain that if T is a triangulation of C and S > 0, then ST is a triangulation 
of 6C. In this way we get the (Ji^i-triangulation of for which mesh{5Ki) = S^/N for 5 > 0. 
Let T be a triangulation of R^ x [0, 1] with the restriction C R^ x {0, 1}, i.e., the vertices 
only lie in R^ x {0, 1}. Then we call this triangulation two-layered. Let Ki denote the 
restriction of the ii'i -triangulation of M'^+-'^ to x [0, 1]. Then Ki is two-layered. Further 
define the {N -\- 1) x {N -\- l)-matrix P = [Sei, . . . , Scn, gn+i] for given 6 > 0. Define 

K,{6) = {{Py\ . . . ,Py^+2) : {y\ . . . G K,}, 

then -^i(^) is a two-layered triangulation of M.^ x [0, 1]. 
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Figure 1 : Illustration of the i^i-triangulation for A'' = 2 
3.2 Some facts about homotopy algorithms 

In this section we want to provide the basic principles of homotopy algorithms. 

Definition 3.3 Let f,g:C—^Dbe two continuous mappings from the topological space 
C to the topological space D. We call f,g homotopic, if there exists a continuous mapping 
1? : C X [0, 1] ^ D, {s,t) ^ ^s,t) with i?(s,0) = f{s) and i?(s, 1) = g{s) for all s e C. We 
call 1? the homotopy from f to g. 

Let C be a nonempty, compact and convex subset of M''^ and assume that / : C — )■ C is 
continuous. Then it follows by Kakutani's fixed point theorem (cf. [1, p. 174]) that there 
exists at least one fixed point of /. To determine any fixed point we use the idea of the 
classical homotopy. Define the continuous mapping /j : C — > C by 

ftis):=il-t)so + tfis), te[0,l] 

with So G C. Then by a further application of Kakutani's fixed point theorem, there exists 
a fixed point of ft for every t G [0, 1]. We start with the constant mapping fo{s) = sq and 
its fixed point so. Assume that — )■ 1 for k — )■ oo, then the sequence of functions (/tj.(-))fceN 
converges even uniformly to /i(-) = /(•)■ Now one can show that the cluster points of the set 
of fixed points st,, of ft,, are just the fixed points of /. Note that in this approach we have to 
extend the dimension of this problem, i.e., we now work in the space C x [0, 1]. 
The numerical procedure for nonempty, compact and convex C C is the following. We 
decompose the space C x [0, 1] in simplices using a suitable triangulation T. Under certain 
conditions there exists a path in this triangulation from an A^-simplex G C x {0} to an 
AT-simplex r* G C x {1} which yields an approximate fixed point of the function /. 
The algorithm that we use here, denoted by SFP-algorithm {simplicial fixed point algorithm) 
for short, follows the path by using the lexicographic pivoting rule from linear programming. 
The advantage is that the so-called degeneration problem (i.e., the path ends up in a circuit) 
cannot occur. We don't want to enlarge on that fact and will only give the definition of 
lexicographically positive matrices. For a detailed description we refer to [31, Chapters 2&;3]. 

Definition 3.4 A row vector is called lexicographically positive, if its first nonzero entry is 
positive. A matrix W is called lexicographically positive denoted by W y if every row 
vector is lexicographically positive. 
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3.3 The SFP-algorithm 

To compute a fixed point of a continuous function cf) : — >■ the SFP-algorithm uses a 
suitable homotopy "9 : M.'^ x [0, 1] — t- and a pivoting method to get from an A^-simplex 
C X {0} to an A^-simplex C x {1} which yields an approximate fixed point of (p. 
For this purpose we have to triangulate the set x [0, 1] suitably. 

Let r be a triangulation of x [0, 1] with the restriction 7^ C x {0, 1}, i.e., T is two- 
layered. We denote elements of x [0, 1] by y = (f i, . . . , VN,t) with v G and t G [0, 1] 
and define the projection onto the first factor pi : x [0, 1] — t- R^, pi{v,t) = v. Suppose 
that the iV-simplex t = {y^ , . . . , y^^^) G . We define the diameter of the projection of 
r by 

diamp(r) :=max{||pi(y^)-pi(y^)|| : i, j e {1, . . . , N + 1}}. 
Moreover, the mesh size of the projection of T is defined by 

meshp(T) := sup{diamp(r) : r G T^}. 

If r = G and r C x {i}, i G {0,1}, then Tp := (pi(y^), . . . is 

an iV-simplex in M^. The collection of all such simplices Tp is denoted by 7^. 
We choose an arbitrary point (c, 0) G x [0, 1] such that (c, 0) lies in the interior of an 
A/"-simplex r° G Tq. Consider the following homotopy mapping -d : x [0, 1] — > M"^ defined 

by 

{}{v,t) = {l-t)c + t(l){v). (15) 

A point y is called a fixed point of -d, if pi{y) = '&{y)- Clearly, (c, 0) is the only fixed point 
of "d in X {0} and any fixed point y of in MJ^ x {1} projects to a fixed point of cf), i.e., 
Pi{y) = (t^{Pi{y))- The concept of labelings establishes a way of studying the relation of the 
triangulation with approximate fixed points of (f). 

Definition 3.5 Let T be a two-layered triangulation o/ x [0,1]. Then we define the la- 
beling rule / : X [0, 1] by 

i{y) = m-pi{y)- (16) 

Let the N -simplex r = (y^ . . . , C be given. Then we call the {N + I) x {N + 1) 

matrix 

''^''^'■={i{y') ::: kA^)) ^^'^ 

the labeling matrix of t. 

The iV-simplex r is called complete, if the system 

L{t)W = Ln+i, H^^O (18) 

has a solution W* G M(^+^)^(^+^). Complete simplices play an important role in the following 
since a complete A^-simplex r C ^ {1} contains an approximate fixed point of (f). In 
addition, by choosing the mesh size of the triangulation small enough we can claim any 
accuracy of the approximate fixed point. 
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Proposition 3.6 Let D he compact and (f) : D C — )■ be continuous. For £>OletS>0 
be such that for all x,y £ D we have the implication \\x — y\\ < 6 ^ \\(l>{x) — (j){y)\\ < e. Let T be 
a two-layered triangulation o/M^ x [0, 1] with mesh(T) < 6 and r = {y^, . . . , y^~^^) C x {1} 
a complete simplex in T with y^ = {v^ ,tj) for all j = 1, . . . , N + 1. Let A G be the solution 
of the system 

L{t)\ = ei, AG M^+^ (19) 
Then v* := Ylf=i^ ^j'"^ approximate fixed point of (f), i.e., \\(f){v*) — v*\\ < e. 

Proof. Since r C x {1} we have tj = 1 for all j = 1, . . . , + 1 and so l{y^) = 
'd{y^) — piiy-') = 4>{v^) — by (15) and (16). Thus (19) is equivalent to 

Af+l N+1 N+1 

(i) ^Xj = l and (ii) Xj(p{v^) = Yl ^i^^ ■ (^0) 

By (20)(i) V* is a convex combination of the v^,. . . ,v^^^, i.e., v* G r. But then we have 
— v-'ll < (5 for all j = 1, . . . , iV + 1 and by continuity of (j) we have ||<;!i('y*) — ^{v^)\\ < £ 
for all j = 1, . . . , A/" + 1. Together this yields 



(20)(n) 



iV^ I 



Af+l A'+l 

mv*)-v*\\ ||(J;A,)</>(^*)-^A,. 

JV+l N+1 

||^A,-</,(t;*)-5;A,>(^;^)|| 
j=i j=i 

. (20)0) 



□ 



To obtain a complete simplex in M.^ x {1} we first characterize the complete simplices. To 
this end we define the graph G(y, E) of all complete simplices as follows. An (A^ + l)-simplcx 
?7 of T is a node, if it has at least one complete facet r. Two nodes are adjacent and connected 
by an edge, if they share a common complete facet. The degree of a node rj is the number of 
nodes adjacent to 77, denoted by deg(77). 

Recall that is the A^-simplex lying on x {0} and containing (c, 0) in its interior. Let 
rf be the unique {N + l)-simplex of T having as its facet. Then we have (cf. [33, Lemma 
4.3.3, Lemma 4.3.4 and Theorem 4.3.5]). 

Lemma 3.7 The N -simplex is the only complete simplex on M.^ x {0}. 

Lemma 3.8 Given the graph G{V,E) defined as above, for each node rj of G{V, E), we have 
(i) if rj has a complete facet lying on x {0} or x {1}, then deg(ry) = 1; 
(a) in all other cases, deg(?7) = 2. 

Theorem 3.9 For the graph G{V,E) defined as above, each connected component ofG{V,E) 
has one of the following five forms 
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(i) a simple circuit (i.e. a path (eo,i,ei^2, • • ■^k-i,k)> A; G N with eo,i = e-k-i^k O'^d ej^j+i 7^ 
ejj+i fori^j andi,j e {l,...,k- 1}); 

(a) a finite simple path (i.e. a path without circuits) whose two end nodes all have a complete 
facet lying on x {1}; 

(Hi) an infinite simple path starting with an {N + 1)- simplex which has a complete facet lying 
on X {1}; 

(iv) a finite simple path which starts with the {N + 1) -simplex rp and ends with another 
{N + I) -simplex having a complete facet on x {1}; 

(v) an infinite simple path starting with the {N + l)-simplex rp . 

From the point of view of computation we are interested in case {iv). In this case we can 
algorithmically go from rp to a simplex rj* G x {1} containing an approximate fixed point 
of (f) by Proposition 3.6. A schematic description is given in Figure 2. 

The Simplicial Fixed Point Algorithm 

Step (0) Set meshp(T) < 5. Let be the unique A'^-simplex of containing 
(c, 0) in its interior. Let rp be the unique {N + l)-simplex in T which has r° as its 
facet. Let be the vertex of rp that is not a vertex of r°. Set A; = 0. 

Step (1) Compute W = L-^{t^) with L from (17). Let Wi denote the i*'* row of 
W. Compute l{y~^) and let q = (1, l'^ {y'^))'^ . Let p = {pi, . . . ,pn+i)'^ = Wq denote 
the coefficient vector of the linear combination l{y~^) = '^^Ji^ pJiy'')- Compute 
C G {1, . . . , A'' + 1} so that the quotient 

^=min<^^ : ph>0, h = 1, . . . , N 1\ 
PC ^ [ Ph J 

is lexicographically positive minimal. Note that C is unique (cf. [33, Theorem 4.2.7]). 
Let r*^"*"^ be the facet of r/*^ opposite y^. If t^~^^ lies on M.^ x {1}, this facet yields 
an approximate fixed point of cf) and stop. If t'''^^ does not lie on x {1}, go to 
Step (2). 

Step (2) Find a simplex rj'^'^^ sharing the facet t^'^^ with rj'^ (which is unique by 
Lemma 3.8), and let y~^ be the vertex r)'''^^ not being a vertex of r'^'^^. Set k = k-\-l 
and return to Step (1). 



Figure 2: The SFP-algorithm 

RemEirk 3.10 In order to guarantee case (iv) in Theorem 3.9 Merrill (cf. [33]) gave a con- 
traction condition that is sufficient for the convergence of the algorithm for a u.s.c. point-to-set 

mapping (j) : — > R^, in particular for a continuous single-valued function cj) : IS.^ — t- M'^, 
see [33, Theorem 4-3.6]. The method of proof is to show that there is only a compact subset 
D C yielding complete simplices, so the path must he finite. Note that this condition is not 
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sufficient for convergence, if we choose cj) : M.^ . In particular the function (f) defined in 

(21) satisfies Merrill's condition (cf. [11, Satz 4-^8]) but we have to impose other conditions 
to guarantee convergence. 



3.4 Using the SFP-algorithm for computing decay points 

Now we want to use the SFP-algorithm to compute a decay point w G J^(rjii) of the monotone 

operator T^-.R^ ^ which satisfies r^(0) = 0. 

In the following the aim is to find a suitable function <p whose fixed points w = (j){w) correspond 
to decay points w G r2(r^), and to show that the SFP-algorithm converges for this choice of 
^. Since r^(0) = and this point yields no information, we have to exclude from being a 
fixed point of 4>. Also, in order to show that complete simplices can only lie in a compact 

subset of M:'y it is desirable to have cj) small for large v. 



Consider the function 



4>{v) = V^{v) l + min<^0, 



defined by 

Kr — 2||t;| 



|f II + Ko 



-|- max{0, Kft — 211^11} e. 



(21) 



Here let > 0, Kr > > and e := X^^^ the A''-dimensional vector of ones. We illustrate 
the components of (j) in Figure 3. 




Figure 3: The components of the function 

Some properties of cf) are as follows: 

(i) is continuous on 

(ii) For large v G 
(hi) It holds (/>(0) 



^ since is continuous on 



it holds (/>(u) < 0. 



KhC ^ 0, i.e., the origin cannot be a fixed point of (j). 

In Figure 4 we illustrate the definition of (p on the positive orthant. To this end we partition 
the positive orthant in five regions: 



I 


= G I 


R^: 


\\v\ 


e [0,Kh/2)}, 


I' 


= lx[0,l]. 


II 


= {v^\ 


: 


\\v\ 


G [Kh/2,Kr/2)}, 


II' 


= nx[o, 1], 


III 


= {ve\ 


: 


\\v\ 


G [Kr/2,Kr + Ko)} , 


III' 


= nix [0,1], 


IV 


= {ve\ 


: 


\\v\ 


G [kt + Ko, lir + kq + S]} , 


IV' 


= IV X [0,1], 


V 


= {vel 


R^: 


\\v\ 


G (Kr + Ko + 5,oo)} , 


v 


= Vx[0,l]. 
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Figure 4: The definition of illustrated on the positive orthant 

The next proposition indicates the relation between fixed points of ^ and decay points of F^. 

Proposition 3.11 Let <^ : ^ 6e defined as in (21) and assume that F^ : ^ 
is monotone and satisfies the small gain condition (14). Let s G M;^ be a fixed point of the 
function (p, i.e., s = ^{s). Then s lies in the set of decay of the function T^, i.e., s G J7(F^). 
Moreover, s £ I. 



Proof. We distinguish between the following two cases for s G Mlj: : 

(i) < < In this case s = (j){s) = F^(s) + (k/j - 2||s||)e » F^(s). It follows that 
s S> F^(s), so s G J^(F^). In particular, s G I. 

(ii) ^ < \\s\\: We have s = 0(s) = F^(s)(l + min{0, ^J^|||W}) < F^(s) but this is a 
contradiction to the small gain condition (14), so this case cannot occur. 

□ 

In the following we will always use the i^i((5)-triangulation. This triangulation has the essen- 
tial advantage that the vertices of an iV-simplex t = {y^, . . . , y^^^) are in the order of M^"'"''^, 
i.e., it holds < . . . < Note that y = {v,t) G M+ x {0, 1}. 

Again, the SFP-algorithm starts with the {N + l)-simplex rp which has the A^-simplex 
G MJ^ X {0} as a facet containing (c, 0) in its interior, where c determines the homotopy 
mapping ■& in (15). Here we choose c G I U II and any approximate fixed point d will also lie in 
lull, see Theorem 3.14. Then the algorithm follows the path of complete A^-simplices. If we 
can show that this path is finite and inside of the positive orthant, then we get, by Theorem 
3.9, that the SFP-algorithm ends up with a (A'^ -|- l)-simplex containing a complete facet on 
X {1}. Proposition 3.6 now tells us that this simplex contains an approximate fixed point 
of (t>. 
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A first rough estimation where the path of complete simpHces can run is given in the next 
proposition. 

Proposition 3.12 Let cf) : he defined as in (21) and assume that : — >■ 

is monotone. Assume that the constant c used in (15) satisfies c G JU // and let r = 
(y^, . . . , y^~^^) be an N -simplex in Y . Then r is not complete. 

Proof. We prove this by contradiction. Assume that r is complete. Then the linear system 

L{t)W = In+i, W^^O (22) 

with L defined as in (17), has a lexicographically positive solution W. We have = {v^,tj) G 
X {0, 1} and by r C V we have ||?;-'|| > Kr + i^o + S, i.e., (t){v^) < for all j = 1,.. .,N+1. 
So we have the following two cases for 

l{y^) = d{v^,tj) - = (1 - tj)c + ti4){v^) - vK 

(i) litj = 1, then l{y^) < 0; 

(ii) If tj = 0, then l{y^) = c - . Since c G I U II and G V for all j = 1, . . . , + 1 there 
exists a component i* G {1, . . . , N} with Cj* < vj* < u^* for all j = 2, . . . , A + 1. 

Together it follows 

Ky^)i* < for all i = 1, . . . , A^ + 1. (23) 

Let Li denote the l^^ row of L and let W"'' denote the m*'^ column of W. Then we have 
G M^^\{0} since W is lexicographically positive. By (23) we have <^ 0. But then 

< in contradiction to (22). So r cannot be complete. □ 

Note that this does not show that the path starting in r)'^ is inside the positive orthant. To 
prove this we have to look at the boundary of the positive orthant. Here we need some 
additional assumptions. Note that for T^, F G (/Coo U {0})''^^''^ is the underlying gain matrix 
(cf. Section 1.3), and by Remark 1.2 the operator is monotone. 

Theorem 3.13 Let (p:R^ be defined as in (21) and assume that the underlying gain 

matrix T is irreducible. Let r = {y^, . . . , y^~^^) be an N -simplex on the boundary of the positive 
orthant. If = ||pi(2/'^^^)|| < + Kr then r is not complete. 

Proof. If T is an A-simplex on the boundary of the positive orthant then there exists an 
index i* G {1, . . . , A} with vj, = for all j = 1, . . . , A + 1. We prove by contradiction that 
r cannot be complete, if < kq + Kr- So assume 

L{t)W = In+1, WyO (24) 

has the solution W* and let A G M^+^ be the first column of W*. Then it follows by (24) 
and using (17) that 

N+l N+1 

^ = 0, Yl = 1- (25) 

j=i j=i 
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The case < ^ yields, using (16), l{y^)i* = (1 - tj)ci* + tj(f>{v^)i* > for all 

j G {1, . . . , iV + 1}, so J2f=i '^jKy^)i* > since A 7^ and c S> (since c lies in the interior 
of a simplex r). But this is a contradiction to (25). 

Now assume < Then it holds l{y^)i* > 0. In particular, l{y^)i* = (1 — tj)ci* + 

tjV^{v^)i*+tj{Kh-2\\v^\\) > for \\v^\\ < ^. Let r G {1,...,A^ + 1} with ti = ... = tr = 
and tr+i = ... = tN+i = 1 as well as p G {1, . . . , AT + 1} with \\vP\\ < ^ and \\vP+'^\\ > ^. 
Then equation (25) is equivalent to 

r AT+l p N+1 N+1 

^ Xjc + ^ A,r^(^^) + ^ XjiKH - 2\\v^\\)e = ^ Xjv^, ^ A,- = 1. (26) 

j=l j=r+l j=r+l j=l j=l 

Since vf, = for all j = 1, . . . , + 1 it follows Aj = for j = 1, . . . , r with f := max{r, p} 
and equation (26) is equivalent to 

N+1 N+1 N+1 

J2 ^j^.iv') = E E = 1- (27) 

j=f+l j=f+l j=r+l 

Now there exists a largest index f G {1, . . . , + 1} with r^(u')i* = and r^(w'')i* > 0. 
Since the are ordered by the ^i((5)-triangulation, it follows by monotonicity of F^, 

T^iv') <...< T,i,/) <...< T^{v^+'), (28) 

and thus Af +1 = . . . = A n+i = which leads to 

f f f 

j=f+i j=f+i j=f+i 

Without loss of generality assume = [v^, . . . ,vl , 0, . . . , 0]"^ with v^j > for J = 1, . . . , Z, 
I < N. Equation (29) implies r^(i'0 = • • • , *i, 0, . . . , 0]"^ with *j > 0, j = 1, . . . , I. But 
then r is of the form 



Til 




(30) 



with Til G M.^^^ and G M^+^ Ox(iV+i 0_ ^his means that T is reducible, a contradiction 
to the assumption. So this case cannot occur. 

Now assume \\v^+^\\ < kq + Kp. Define t{v^) := T^{v^) (l + ttS^) if Ib^ll > ^- 



-MV- 7 -r y II ^ 2 

||^^"'|| < /«o + it holds T{v^)i^ > and Tn(v^)k = 0, if and only if r{v^)i^ = 0. The same 
argumentation as above provides 

E mv') = E E = 1- (31) 

j=r+l J=r+1 j=r+l 

With = K, . . . , t;f, 0, . . . , 0]T, > for j = 1, . . . , / it follows f(t>0 = 0, 0]^ 

with ij > 0, j = 1,..., I. All in all we get t{v^) < Tu,{v^) = [*i, . . . , 0, . . . , 0]"^ with 
*j > *j > 0, j = 1, . . . , Z. But then T is of the form (30), a contradiction to the assumption. 

□ 
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In other words Theorem 3.13 provides that no A^-simplex r G I' U II' U III' lying on the 
boundary of the positive orthant can be complete. So it remains to show that the path 
starting in rf cannot enter the set IV'. For this purpose we show in the next theorem that 
the path of complete simplices runs inside of the region which is painted dark grey in Figure 
5. To prove this we demand an upper bound for the feasible size of 5. 
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Figure 5: -fCi ((5)-triangulation and the maximum region of the path (r-simplices are 1- 
dimensional and ry-simplices are 2-dimensional) 



Theorem 3.14 Let (f) : 'R'^ — >■ M-'^ he defined as in (21) and assume that : 
satisfies the small gain condition (14). Furthermore assume that the underlying gain matrix F 
does not contain any zero row, i.e., T^{e)i ^ for all i = I, . . . , N . Then there exists a 5 > 
such that for all simplices t = {y^,..., y^~^^) C /// with := y^ + [6, . . . , S, 1]^ G iV it 

follows that T cannot he complete. 

In particular, any approximate solution d satisfies d G /U II. 

Proof. Simplices satisfying r C III' with y^^^ := + [(5, . . . , b, V\' € IV' are marked as 
black dotted lines in Figure 5. We show that such a simplex r cannot be complete, i.e., the 
system 

L(t)VF = 7jv+i, VF^O (32) 



has no solution. First it holds for all s G with ||s|| < Kr + ko that 
r^(s) < max{r^(i;) : t; G A ||w|| = Kr + Ko} 

< r^([Atr + Ato, • • • , /«r + /«o]^) =: r;; 

Choose A; G N such that 



^max 



^ -V'T''<C (33) 



2A;-1 
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and set (5 > such that 

6 < mm < > . (34) 

I 2ViV 2kVN J 



Now any approximate solution c' can only be in I U II and not in III since 6 < '^^^/^ ' This 
follows by Proposition 3.11 because any fixed point s* of has the property ||s*|| < ^ and 
then 

II 'II / II *ii , I. (i? ixw II *ii I ^ /^^/^ '^'1 _L i^T 
c < \\s + mesh„(Ai(())) = s + ov N < — H -j^^v N = —. 

^ 2 2Vn 2 

The idea now is the following: Let r G {1, . . . , N} with = and tr+i = 1. Then it holds 



L(r) 



1 ... 1 1 ... 1 



Moreover, G III. So there exists at least one index I G {1, . . . , + 1} with (c — v^)i < 0. 
Then 

{c-v>)l < for j = l,...,r. (35) 

Now (f){v) converges to zero, if v tends to the boundary of IV, i.e., iiv v with \\v\\ = nr+i^o- 
So the aim is to get c/){v) as small as {^{v^) — v^)i < holds for all j = r + 1, . . . , iV + 1. 
Note that the function 

5 : M+ ^ M, a >-> 1 + — 

a + kq 

is strictly decreasing for k,q > 0, Kr > 0. Prom the relation 1;^+^ = + [6, . . . 5]^ we get 

+ < 11^^^"^^ II = 11^;^ + [(5, . . . , 5^ \\ < \\v^ \\ + 6^/N. 

Under this assumption it follows for all a > ky + kq — S-\/N 

, - 2{Kr + K0- 6VN) 
g{a) < 1 + 



Kr + — SV N + Ko 
Kr + 2ko - S^/N + Kr- 2Kr - 2ko + 25\/N 



KT + 2kq - 5\/N 



SVN (34) SVN 1 
< 



Kr + 2ko - SVN 2k6VN - 5VN 2k -I 
Together with (33) it follows for j = r + 1 , . . . , + 1 

V ||i;-?|| + Ko J \ + f^o J 

1 . (33) 

Altogether with (35) it follows for all j = 1, . . . , iV + 1 

Wl = ((1 - tj)c + tj4>{v^))i - vi <ci- vi < 0. 

Let Li denote the Z*'* row of L and let denote the m*^ column of W. Prom the above 
consideration it follows {Li)'^ <g and from LiW^ = 1 we get > 0. But then it follows 
LiW^ < 0, a contradiction to LiW^ = according to equation (32). So r is not complete. 

□ 

Now we can deduce the following main theorem. 
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Theorem 3.15 Let (p be defined as in (21) and assume that T is irreducible and that the 
operator : M.^ — )■ M^, deduced from the gain matrix F, satisfies the small gain condition 
(14). Let 6 > be chosen as in (34) with k & N according to (33). Then the simple path 
starting with rp is finite and the SFP-algorithm converges to a decay point s G S7(r^). 

Proof. The dark grey painted region in Figure 5 is compact. The path of complete simphces 
starts in the interior of this region. Theorem 3.13 and Theorem 3.14 now show that under the 
above assumptions the path starting with rj^ cannot leave this region. So the path remains 
in this region. Since the region is compact there exist only finitely many simplices in this 
region and we arc in the situation of Theorem 3.9 (iv) . So the path is finite and ends up 
in a simplex r G M.^ ^ {1} which contains an approximate fixed point of (p by Proposition 
3.6. So refining the triangulation leads to the convergence of the SFP-algorithm to a fixed 
point s of 4>. Since satisfies the small gain condition (14) it follows by Proposition 3.11 
that the fixed point s of (p lies in the set of decay Q{Ti^). So the SFP-algorithm converges 
to a decay point s € ri(r^). □ 

3.5 Improvement of the algorithm 

To summarize implementation details we give some suggestions for the choice of c, S, and the 
constants kq, kh and Kr for a given function F^ of dimension N. 

Suggestions for the choice of Kh,Kr,Ko 

Theorem 3.11 says that a fixed point s = (f){s) can only lie in region I. Since (f){s) = F^(s) -|- 
(k/i — 2||s||)e for s G I we may expect the fixed point s to have a norm near Kh/2. So we 
choose Kh as the double size of the norm of the desired fixed point. 

Several computational experiments have shown that values for /cp near Kh and kq small will 
probably lead to small computing times. So we give the suggestions Kr = k/i + 1 and kq = 1- 
Note that for smaller values the regions II' and III' are small and so the path tends to leave the 
region I' U II' U III' more often. This leads to more pivoting steps and so to longer computing 
times. 

Suggestion for the choice of c 

If we have no advance information about the location of the fixed point we choose c by default 
as c = 0-99^^ [1 . . . 1]^. The norm of c then is ||c|| = 0.99k/i/2, thus near k/j/2 where we 
expect the fixed point. In addition no direction is preferred. 

In some cases we have some information about the approximate location of the decay point. 
Then we can use this information by using this expected point as c (if it lies in I U II) to arrive 
smaller computing times. 

Suggestion for S and the refinement of S 

The choice of 5 > as in (34) is one that leads to provable convergence but we have seen in 
experiments that this choice leads to longer computing times. So we will ignore the choice 
of (5 as in (34) and give another suggestion. To ensure that the algorithm converges stop the 
iteration, if the path leaves the region I' U II' U III', and start again with the same starting 
point and a refined, i.e., smaller 6. Since the dimension N of the operator F^ can get large it 
is advisable not to choose S too small. Since the simplex has a diameter of \/iV S we also have 
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to include the dimension N into the choice of 6 such that the path doesn't leave the region 
I' U II' U III'. Our suggestion therefore is 5 = 

The algorithm refines d until the desired accuracy is reached. So it is important that we 
do not only choose S suitably, but even determine the refining sequence {(5fc}fegN such that 
the approximation is quite good. Saigal [25, Section 5] gave such a sequence and showed 
that the algorithm converges quadratically, if we assume in addition that (f) is continuously 
differentiable and its derivative is Lipschitz continuous. 

In our case we are content with the refining factor ^, i.e., 5k+i = The simple reason is 
that the algorithm stops, if it finds a point in the set of decay. Again numerical experiments 
suggest that we do not have to refine often to find a decay point. 

SumniEiry 

We want to summarize our suggestions. Let nomi> denote the desired norm of the decay 
point. Then we have the suggested values to be computed as 

Kh = 2nomi, Kr = + 1, kq = 1, c = 0.99 — 7=e, (5 = — ^. 



Remark 3.16 Note that if the SFP-algorithm does not converge for 6 as in (34), then the 
small gain condition ^ id cannot hold on the whole region JU IIU IIIU IV. So we have to 
choose a smaller norm> 0. 

4 Examples 

In this section we give two examples. First note that in [24] an algorithm is developed to 
compute decay points, which is derived from a homotopy algorithm due to Eaves [9]. In [23] 
an example is given and decay points are computed with the algorithm from [24], referred to 
as Eaves algorithm. We state the principle ideas of this article, pick up the results given by 
Eaves algorithm, and compare them to those of our SFP-algorithm. 

The second example concerns about a biochemical control circuit model which leads to a non- 
linear gain matrix F. We give a general example of monod kinetics and state some conclusions 
about the input-to-state stability of this control circuit model. At the end we consider a per- 
turbed system and use the methods presented here to check the local input-to-state stability 
numerically. 

4.1 Quasi-monotone systems 

The motivation for this example was the article of Riiffer et al. [23]. Therein a nonlinear 
system is given and decay points are computed with Eaves algorithm from [24]. For this 
purpose a nonnegative matrix P G M.^^^ with spectral radius p{P) < 1 is constructed for 
given dimension N. By Perron-Frobenius theory it follows that the matrix A := —In + P 
then has spectral abscissa a{A) := max{ReA : A is an eigenvalue of A} = —1 + p{P) < 0. So 
the matrix A is Hurwitz with negative diagonal entries and nonnegative off-diagonal entries. 
Now we define a smooth coordinate transformation S : — >■ M.^ by 



S{v)i = 



< 



V. 




if > 1 

if Vi e [-1, 1] 

if Vi < -1 
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It holds 5(0) = and S{R^) = R^. The mapping S : ^ is a monotone operator. 
Then the systems 

v = S'{S-\v))AS-\v)=:g{v) (36) 

and 

i = Az=: h{z) (37) 

are equivalent under a nonlinear change of coordinates. Let v* be any decay point for the 
function g in equation (36). With it z* := S~^(v*) is a decay point for the function h in 
equation (37). Wc want to pick up the associated run times and numbers of iterations and 
compare them with those of the SFP-algorithm. 

The following results correspond to matrices P G M^^^ with positive entries in [0, 1] generated 
by a numerical approximation of the uniform distribution, and 30% of those are set to zero. 
Then a{A) = —0.2. The numbers are averages over 100 simulations. Here we assumed 
norm = 10, i.e., the norm of the desired decay point v* is \\v*\\ ~ 10. In Table 1 the results 
of [23] are listed. In Table 2 we give the results of the SFP-algorithm. In addition, we tested 
the SFP-algorithm even for large N. 



N 


run time 


# iterations 


5 


0.11465s 


267.62 


10 


0.64855s 


2059.65 


15 


1.7833s 


5505.78 


25 


7.987s 


19742.84 



Table 1: Results of Eaves (Kl)-algorithm from [23] for norm = 10 



N 


run time 


# iterations 


simulations 


5 


0.0277s 


20.9 


100 


10 


0.0415s 


34.5 


100 


15 


0.0618s 


72.3 


100 


25 


0.1710s 


187.8 


100 


50 


1.180s 


688.4 


100 


100 


13.22s 


2711.9 


50 


150 


78.35s 


6614.3 


10 


200 


273.6s 


11243.8 


10 



Table 2: Results of the SFP-algorithm for norm = 10 



Note that the run times and iterations can only be compared relatively since the simulations 
are executed on different computers. Nevertheless, our run times are considerably lower and 
even for relatively large dimensions we are able to compute decay points in a quite acceptable 
run time. 

In Table 3 we give run times and iteration numbers for norm = 1000. One can see that we 
have a relatively small increase of iteration steps and therefore of run times despite a quite 
larger norm. This is a consequence of the fact that we choose the size of S, and with it the 
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N 


run times 


7^ iterations 


simulations 


5 


0.0451s 


61.6 


10 


10 


0.0680s 


62.5 


10 


15 


0.0879s 


106.6 


10 


25 


0.2977s 


317.0 


10 


50 


2.141s 


1097.1 


10 


100 


25.07s 


3991.3 


10 


150 


253.2s 


9214.9 


10 


200 


542.0s 


16252.1 


10 



Table 3: Results of the SFP-algorithm for norm = 1000 



mesh size of the starting triangulation in dependency of the norm norm [6 = ^°^™ ). 

That is why the algorithm gets close to the desired decay point in few steps. 

4.2 A biochemical control circuit model 

We consider the following biochemical control circuit model similar to [27] , 

xi{t) = g{xN{t)) - aixi{t) + u{t) 

Xi{t) = Xi-i{t) - aiXi{t), i = 2,...,N, (38) 
x{t) = [xi{t),...,XNit)V eR^, 

with ai> constant for alH = 1, . . . , A'", it G L°°{[0, oo); M) and g : M_|_ M_|_ a continuously 
differentiable function with g{x) > for all x > 0. In contrast to [27] we added an external 
input u and do not assume g to be bounded, but we demand the following assumption, which 
was introduced in [17]. 

Assumption 4.1 There exist x*j^ > 0, > and A G (0,1) with ax*j^ = g{xy) and a = 
YljLi dj such that 



K + x*^ 
K + x 

Remark 4.2 The function 



X < a ^g{x) < x%- + X\x — for all x > 0. (39) 



bx 

g{x) = ^-p^, x>0, 6, c> (40) 

models the growth rate of cells or micro-organism and is further known as monod kinetics. Let 
a > be arbitrary and b,c> such that b > ac. Then assumption 4.1 is satisfied. This can 
easily be seen by setting x*^ := > 0, iiT = c > and A := . 

Following similar calculations as in [17] we get the following result by setting the gains as 
li,j{s) = d forjVA^and 7i,7v(s) := ^ (in (l + ^(exp(V2^) - 1)))', 

li,j{s) = foriT^i-l, 7i,i_i(s) := ^(ln(l + C(exp(V2i)-l)))', 

with e G (max{^^,A},l) and C e (1,6"^). 
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Theorem 4.3 [11, Satz 5.5] Consider the system (38) with g defined as in (40) with b,c> 0. 
If b > ac with a = Y[f=i > then the equilibrium solution x* := [x\, . . . ^ with 

x*^ := ^ and x* := (Uf^i^ x*^ for i = 1, . . . , N - 1 is ISS on M^\{0}. 

4.2.1 A perturbed biochemical control circuit model 

Theorem 4.3 is a nice theoretical resuh. But in apphcations we are always faced with perturbed 
systems. In this section we want to check the local input-to-state stability of a perturbed 
system with the methods developed in this work. 

For this purpose consider the graph in Figure 6. The black arcs describe the real couplings 
of the biochemical control circuit model and the grey arcs describe the perturbations. The 
underlying system to this coupling graph is given by 

xi{t) = g{xs{t)) - aixi(t) + u{t)+eu{xi{t)) 

= Xl{t) - a2X2{t)+i22{x2{t)) (41) 

xsit) = X2{t) - a2,xz{t)+eiz{xz{t)) + ezi{xi{t)). 




Figure 6: The coupling graph of a perturbed control circuit model 

For this system the functions e describe the perturbations. Further let oi = 2, as = 1, 03 = 3 
and a, = 04 • 02 • 03 = 6, and g : — t- R+ be the monod function given in (40) with 6 = 8 and 
c = 1. Then the associated gain matrix T is of the form 



711 
721 722 
731 732 733 



713 





Here let 



7i3(s) := ^ (ln(l + 0(exp(V2^) 
721 (s) := ^ (ln(l + C(exp(V2^) 
732(s) := ^ (ln(l + C(exp(x/2^) 



I)))- 



be the gain functions from the previous subsection with 

K 



(42) 
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By Remark 4.2 it follows K = 1, X = ^ and = 1/3, i.e., 9 G (|, !)■ Here we assume 
9 = 0.8 and = 1.1, where (42) is satisfied. Since b = 8 > 6 = ac and a = 6 > 0, it follows by 
Theorem 4.3 that the equilibrium solution x* = [ 2 1 1/3 ]^ of the unperturbed system 
(en = £22 = £31 = £33 = respectively 711 = 722 = 731 = 733 = 0) is ISS on R^\{0}. 
The perturbed gain functions of the system are given as 

7ii(s) := 0.001s, 73i(s) := O.OOSs^, 722(5) := 0.001s°-^ 733(5) := O.OOls^. 

For the monotone aggregation function fj, = ^ we get the monotone operator F^. Now we 
apply the SFP- algorithm with the suggestions in section 3.5 and norm = 12 to F^ and get the 
decay point as 





' 6.54 " 




" 6.527 " 


w = 


6.90 


> 


6.886 




7.33 




7.325 



Since limk-^oo^^iw) = 0, Theorem 2.2 is applicable, so the perturbed system (41) is locally 
ISS. 

Finally we illustrate the first and third components of the path a starting in w as defined in 
(11) in Figure 7. Recall that a is obtained as a linear interpolation ( — ) of the points T^{w) 
(•), where we plot this for k = 1, . . . , 1000. Indeed, straightforward calculations show that the 
line from to w ( — ) is not contained in the decay set 17 (F^). Although the difference between 
the path and the straight line appears to be negligible we see that without the numerical effort 
of computing the ^^{w) no O-path is obtained. For a better view we enlarged one region. 



Function SIGMA 




Figure 7: Components of the path a from to w 
4.2.2 A higher dimensional test 

In this section we consider system (38) and want to study the computational effort for higher 
dimensions. Let G N denote the dimension of system (38) with = ^i^, i = 1, . . . ,N, and 
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the function g be defined as in (40) with c = 1 and b = 2N, then Assumption 4.1 is satisfied 
with 

N 

a = l[ai = N + l, x% = ^, K = l, A=^. 

1=1 

Now by Theorem 4.3 this system is ISS on M^^\{0}, if 

^G(A,1) = (^,1) and Ce(l,^-^). 

In Table 4 we tested the computational effort for higher dimensions for noriii=12. Note that 

1 

for large TV we have 9 ^-i near 1, so we have only a small range in choosing ^ such that the 
system (38) still is ISS on R:^\{0}. So we guess that the decay set will be very thin and so 
the decay points are harder to reach, resulting in longer computing times. This can be seen in 
Table 4. Note that the run times include checking that the sequence ^|^(^^;) is a zero sequence. 

The counter k_step indicates the first ^ G N such that ||r^(it;)|| < 10~^. 



N 


e 


C 


run times 


7^ iterations 


k_step 


10 


0.75 


1.020 


0.30s 


134 


1215 


50 


0.75 


1.003 


4.99s 


1405 


4501 


70 


0.75 


1.002 


1.72s 


74 


5911 


90 


0.75 


1.002 


57.61s 


8426 


10257 


110 


0.70 


1.002 


105.25s 


9632 


9888 


150 


0.70 


1.001 


532.43s 


22856 


8961 


200 


0.70 


1.001 


2168.18s 


52752 


12656 



Table 4: Results of the SFP-algorithm for norm = 12 



5 Conclusions 

In this paper we have presented a homotopy algorithm, that is suitable for the computation 
of decay points of gain operators which are crucial in checking the local input-to-state stabil- 
ity. The algorithm is proved to converge in a semi-global fashion, provided the mesh size of 
the underlying triangulation is sufficiently small, but experiments suggest that the result is 
conservative and that larger mesh sizes are frequently sufficient. The algorithm improves on 
a previous simplicial algorithm. The advantage of such algorithms is that they can be used to 
analyze networks with quite general small gain formulations whereas other approaches rely on 
special structure like linearity of the gain operator or the use of maximization as the monotone 
aggregation function. In future research we intend to further develop numerical techniques for 
small gain results and explore relevant examples. 
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